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ABSTRACT 

Supermassive black holes are not only common in the present-day galaxies, but billion 
solar masses black holes also powered z > 6 quasars. One efficient way to form such 
black holes is the collapse of a massive primordial gas cloud into a so-called direct 
collapse black hole. The main requirement for this scenario is the presence of large 
accretion rates of > 0.1 Mo/yr to form a supermassive star. It is not yet clear how and 
under what conditions such accretion rates can be obtained. The prime aim of this 
work is to determine the mass accretion rates under non-isothermal collapse conditions. 
We perform high resolution cosmological simulations for three primordial halos of a 
few times 10^ Mq illuminated by an external UV flux, J 21 = 100- 1000. We find that a 
rotationally supported structure of about parsec size is assembled, with an aspect ratio 
between 0.25-1 depending upon the thermodynamical properties. Rotational support, 
however, does not halt collapse, and mass inflow rates of ~ 0.1 Mq/yc can be obtained 
in the presence of even a moderate UV background flux of strength J 21 > 100. To 
assess whether such large accretion rates can be maintained over longer time scales, we 
employed sink particles, confirming the persistence of accretion rates of ~ 0.1 Mo/yr. 
We propose that complete isothermal collapse and molecular hydrogen suppression 
may not always be necessary to form supermassive stars, precursors of black hole 
seeds. Sufficiently high inflow rates can be obtained for UV flux J 21 = 500 - 1000, at 
least for some cases. This value brings the estimate of the abundance of direct collapse 
black hole seeds closer to that high redshift quasars. 

Key words: methods: numerical - cosmology: theory - early Universe - high redshift 
quasars- black holes physics-galaxies: formation 


1 INTRODUCTION 


The observations of high redshift quasars at z > 6 reveal the 
presence of supermassive black holes of about a few billion 
solar masses in the early universe l^netalj_200^ 200^ 
Willott_et_alJ_ 2010l: iMortlock et al.l I 2 OII ; Venemans et al.l 


20131: IWu et al.ll2015ir From a theoretical perspective, the 


assembly of such massive objects within the first billion years 
after the Big Bang is challenging. How and under what con¬ 
ditions they were formed is a question of prime astrophysical 
interest. 

Various models for the form ation of supermassive black 
hole seeds have been propo sed jRees 1984; IVolonterilI 2 OIOI : 


IVolonteri fc Bellovarvll2012l : iHaiman 


2012fl . One potential 


pathway is the collapse of a primordial, i.e., metal-free, “nor¬ 
mal” star into a stellar-mass black hole. Earlier numerical 
simulations performe d to study the f o rmation of the first 
generation of stars jAbel et al.l I 2 OO 2 I : iBromm et al. I I 2 OO 2 I: 
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lYoshida et al.l [ 2 OO 6 II suggested that the masses of the first 
stars are of the order of a few hundred solar. However, re¬ 
cent simulations found that the protostellar dis k forming in 
the minihalo fragments into multiple clumps dciark et al.l 
201l|_ Greif_et_^ 201^ Stacv et al.ll201 iLatif et~ 2013bl : 


Latif fc Schleicherll2015 1 and leads to the formation of mul¬ 
tiple stars. In the case of a single star per halo, they may still 
reach up to a thousand solar masses at z=25 llHirano et al.l 
l2014h . but feedback from the blac k hole itself shuts its own 
accretion and limi ts its growth (Ijohnson fc BrommI 20071: 
lAlyarez et al. li^. Perhaps, as suggested bv iMadau et ^ 
I 2014I1 . they may still grow under prolonged episodes of 
super-Eddington accretion. 

The second possibility could be the collapse of a dense 
stellar cluster into a massive black hole due to relativis¬ 
tic instability (IShapiro fc Kan3ll987ll or the formation of 
a black due to the stellar dynamical processes in the first 
stellar cluster ([ Peyecc hi fc Volonteril l2009l : iGoswami et al.l 
|2014 iKatz et al.l 2015l l or even the core collapse of a dense 
cluster of stellar mass black holes leading to the formation 
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of a m assive black hole seed d Davies et al. et alJ 

l2014h . The expected mass of a seed black hole from these 
scenarios is a few thousand solar masses. 

A third mechanism could be the collapse of a pro- 
togalactic gas cloud into a massiv e central object, a 
so-called direct collapse black hole (ILoeb fc Rasid fl9^: 


Bromm^ Loebl20^:lBegelman et al.l2006| : lHosokawa et al] 

2 OI 2 I : ISchleicher et al.l I 2 OI 3 II . A detailed discussion of 


these mechanisms is given in reviews (IVolonteril l20ld : 
IVolonteri fc BellovarvI I 2 OI 2 I : iHaimanl 1201211 , and we sum¬ 
marize here the main features of relevance for this in¬ 
vestigation. The particular direct collapse scenario we ad¬ 
dress here requires the pr esence of large accretion rates of 


about > 0.01 - 0.1 MQ/yr ( 

BeeelmanI 2010l: Hosokawa et al.l 

l2013l:ISchleicher et alJl201,‘- 

: Ferrara et al.ll2014l). This is be- 


cause such accretion rates are the main pre-requisite to 
rapidly build a supermassive star of about 10“* - 10^ Mq, 
which would eventually colla pse into a black hole retain¬ 
ing most of the star’s m ass dBaumgarte fc Shapiro! 1 19991 : 
IShibata fc Shaoir^ I 2 OO 2 II . These large accretion rates can 
either be obtained through dynamical processes or thermo¬ 
dynamically, by keeping the gas warm (i.e, with a higher 
Jeans mass) to avoid fragmentation and star formation. 
iBegelman et al.l d2006l ) proposed that self-gravitating gas 
in dark matter halos may loose angular momentum via 
‘bars-in-bars’ instabilities and efficiently assemble a central 
massive star. On the other hand, thermodynamical direct 
collapse requires the suppression of molecular hydrogen to 
avoid fragmentation which leads to a monolithic isothermal 
collapse. In both cases, the halo has to transfer angular mo¬ 
mentum efficiently to avoid collapse to b e halted by t he an - 
gular momentum barrier. For instance, IChoi et aP l|20l4 ) 
have found that angular momentum in isothermal conditions 
can be transported by non-axisymmetric perturbations. 

The conditions in primordial halos are ideal for the for¬ 
mation of massive stars due to the lack of efficient coolants 
in the absence of dust and metals. However, the presence of 
molecular hydrogen may lead to efficient fragmentation and 
star formation. UV stellar radiation can quench Ht forma¬ 
tion and consequently fragmentation can be avoided. In fact, 
numerical simulations show that, in principle, in the absence 
of molecular hydrogen, a metal-free halo collapses isother- 
mally via atomic line cooling. Consequently, a supermassive 
star of about 10 ^ Mq co uld be assembled v ia large accretion 
rates of about 1 Mq/vi llLatif et al.l I2013al l3: 1 Johnson et al.l 
1201.31: iRegan et al.ll20l4) . 

The above scenario requires a high UV flux to disso¬ 
ciate molecular hydrogen, which depends on the radiation 
spectra of stellar populations. We define Jem as the value 
of the UV background flux above which full isothermal col¬ 
lapse occurs because molecular hydrogen formation is sup¬ 
pressed. Jcrii depends on the spectrum of the stars contribut¬ 
ing to the background. For hot stars with Trad = 10^ K, such 
as Population III stars (P op III), Jcru is > 10“* (in units of 
J 71 = 10 ~^‘ erg/cmVs/H z/sr lOmukaill 200 d : IShang et al.ll 2 O 10 l : 
lUatif et aI.ll20lA l2015h . The second generation of stars is, 
however, favoured as source of UV radiation for this sce¬ 
nario due to their longer lives and higher abundance. If such 
stars have a softer spectrum, with Trad = 10“* K, the Jent is 
400-700, much reduced compared to the Pop III case. Re¬ 
cently, however, the spectra of Pop II stars were computed 
employing the stellar synthesis code STARBURST, finding 


that realistic stellar spectra of Pop II stars can be mimicked 
by adopting r adiation spectra with temperature between 
10“* - 10^ K lls^ imura et al.l l2014l : lAgarwal fc Khochfarl 
l2015h . lUatif et ahl ( 2 OI 5 I 1 found that the value of UV flux 
necessary to obtain full isothermal collapse for realistic spec¬ 
tra of Pop II stars is a few times 10“* in units of Ui and be¬ 
comes constant for Rad between 2x 10“* - 10^ K. lRegan et al.l 
ll2014l) considered an anisotropic monochromatic radiation 
source and found that a small amount of H 2 forms due to 
self-shielding even under very strong flux. They further pro¬ 
pose that J 21 = 1000 may be sufficient to form a supermas¬ 
sive star. Similarly, high values of J^m are required for com¬ 
plete isot hermal collapse when X-ray ionization he ating is 
included llUatif et al.ll201^ : llnavoshi fc Tanakall20l3 ). 


Within a scenario that requires full isothermal collapse, 
the number density of direct collapse black holes strongly 
depen ds on the value of Jerit- The estimates of lPiikstra et al.l 
(l2014h suggest that 7^™= 1000 seems sufficient to produce 
the observed number of quasars at z > 6. However, the black 
holes density drops by about a few orders of magnitude be¬ 
low their observed abundance for Jcm> 10“*. This makes the 
direct collapse model infeasible by making their sites ex¬ 
tremely rare. Hence, it becomes necessary to explore whether 
high accretion rates can be maint ained under non-is othermal 
conditions. In fact, the work of iLatif et al.l (|2014| ) for soft 
idealised radiation spectra of 10“* K suggests that it might 
be possible to form massive stars even in the presence of a 
moderate amount of H 2 . However, it is still unclear whether 
for realistic Pop H spectra high accretion rates of 0.1 Mo/yr 
can be maintained for fluxes lower than J„./, and what value 
of the UV flux may be sufficient for this purpose. If so, are 
these sites abundant enough to reproduce the number num¬ 
ber density of quasars at z > 6? We address these questions 
in the present work. 


In this study, we explore the formation of supermas¬ 
sive stars under non-isothermal collapse conditions and de¬ 
termine whether the large inflow rates of 0.1 Mg/yr neces¬ 
sary for their assembly can be obtained. We perform three- 
dimensional cosmological simulations to study the collapse 
of massive primordial halos of a few times 10 ^ Mq form¬ 
ing at z >10 and vary the strength of the background UV 
flux emitted by Pop H stars. In addition to this, we employ 
a detailed chemical model which includes all the necessary 
thermodynamical processes and solve them self-consistently 
with the cosmological simulations. We derive the dynamical 
and thermodynamical properties of these halos. The tempo¬ 
ral evolution of the mass inflow rates is computed beyond 
the initial collapse by employing sink particles. Our results 
show that one may not necessarily need complete isothermal 
collapse to form a supermassive star. This work has impor¬ 
tant implications for the observed abundance of quasars at 
z >6 and provides a potential pathway for the formation of 
supermassive black holes at high redshift. 

This article is structured as follows. In section 2, we 
describe the simulations setup and numerical methods. We 
present our main results in section 3. In subsections 3.1, 
3.2 and 3.3, we show the structural properties of halos and 
their time evolution, the expected mass inflow rates and the 
mass estimates of sinks. In section 4, we summarise our main 
findings and conclusions. 
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Table 1. Properties of the simulated halos. 


Model 

No 

Mass 

Mo 

Redshift 

z 

jcnX. 

Dad = 2 X 10“* K 

spin parameter 

A 

A 

5.6 X 10'' 

10.59 

20000 

0.034 

B 

4.06 X 10^ 

13.23 

40000 

0.02 

C 

3.25 X 10^ 

11.13 

50000 

0.03 


2 NUMERICAL METHODS 


The simulations presented in this study ado pt the open 
sourc e three dimensional code enzo version 2.4 j jBrvan et al.l 
I 2 OI 4 II . Enzo is an adaptive mesh refinement (AMR), par¬ 
allel, grid based, cosmological hydrodynamical simulations 
code. It utilises the message passing interface (MPI) library 
for parallelisation and is well suited for the current simula¬ 
tions. The hydrodynamical equations are solved using the fi¬ 
nite difference scheme by employing the piece-wise parabolic 
method (PPM). The dark matter (DM) evolution is handled 
by the particle-mesh technique commonly used in the Eule- 
rian codes. 

Our simulations start with cosmological initial condi¬ 
tions at z=100 and are generated from the ’’inits” package 
available with enzo. The parameters from the WMAP seven 
years data release ar e used to generate the initial conditions 
l|,Iarosik et al.l[^ll^ . The computational volume has a co¬ 
moving size of 1 Mpc/h and periodic boundaries both for 
hydrodynamics and gravity. We first run simulations with 
a uniform grid resolution of 128^ cells (also 128^ DM parti¬ 
cles) and select the most massive halos forming in the box 
at z=15. The simulations are restarted with nested grid ini¬ 
tial conditions, one top and two nested refinement levels 
each with a grid resolution of 128^ cells. We employ in to¬ 
tal 5767168 DM particles to solve the dark matter dynamics 
which results in DM resolution of about 600 Mq. In addition 
to this, we add 18 dynamical refinement levels during the 
course of simulations which provide an effective resolution 
of about 1000 AU in the central 62 kpc comoving region. Our 
refinement criteria is based on the gas density, the particle 
mass resolution and a fixed resolution of 32 cells per Jeans 
length. The cells are flagged for refinement if their density 
is greater than four times the cosmic mean or the particle 
density is 0.0625 times here r is the refinement factor 

(i.e. 2), 7 is the refinement level, and a = -0.3. To avoid spu¬ 
rious numerical affects, DM particles are smoothed at 12th 
refinement level, which corresponds to a resolution of about 
2.7 pc in comoving units. 

This technique of smoothing DM parti cles has been em¬ 
ployed in numerous high resolution studies llWise et al .I2OO8I: 


Turk_et_^ 20]j_ Latif et al.ll2013al . However, we note that 


Regan et al.l ll2015ll have recently found that coarse DM res¬ 


olution leads to spurious numerical effects such as artificial 
heating of the gas. This effect is stronger when no smoothing 
is applied to DM particles, and even smoothing is not suf¬ 
ficient in some cases, e.g., when the DM mass resolution is 


* http://enzo-project.org/, changeset:48de94f882d8 


above 800 Mq for their cases with lower background UV flux. 
They suggest that to obtain numerically converged results 
the baryonic mass within the core should be approximately 
100 times the DM mass. To asses whether the baryonic core 
is properly resolved in our simulations, we show the ratio 
of enclosed DM to baryonic mass in figure [T] for halo A, as 
a representative case. Our results show that the DM mass 
dominates above 100 pc for J 2 i= 500 and 1000 cases while 
for 121 = 100 it extends down to 10 pc. Below these scales, 
the ratio of enclosed masses (MoM/Mba^yon) declines below 
unity for all BUV flu xes, with the s ame b ehaviour as in the 
well resolved cases in iRegan et al.l ll2015ll . This shows that 
at small radii the baryonic core mass is about 50-100 times 
the DM mass in our simulations, it is sufficiently resolved 
and not affected by spurious DM resolution. Similar trend 
is observed for other two halos. 


To solve the thermal evolution of the gas along with 
cosm ological simulation s, we employ the KROME pack¬ 
age (iGrassi etyjJj JgOl^. Our chemical model is described 
in detail in Latif et ahl ll2015h and here we briefly sum¬ 
marise its main features. The rate equations of H, H+, He, 
He+, He++, e“, H“, H 2 , H/ are self-consistently solved with 
the hydrodynamics. We employ a uniform isotropic back¬ 
ground ultraviolet (BUV) flux with T^ad = 2 x 10"^ K emit¬ 
ted by Pop H stars in units of J 21 = 10“^' erg/cm^/s/Hz/sr 
and turn it on at z=30. Our model includes H 2 forma¬ 
tion, H 2 photo-dissociation, H“ photo-detachment, H 2 col- 
lisi onal dissociation and t he Hp self-shielding fitting formula 
of IWolcott-Green et al.l (l201lll . We include various cool¬ 
ing/heating mechanisms such as molecular hydrogen line 
cooling, cooling due to the collisional-induced emission, cool¬ 
ing/heating due to the three-body reactions, chemical cool¬ 
ing/heating and atomic line cooling (i.e. cooling due to colli- 
sional excitatio n, collisi o nal ion isation, radiative recombina¬ 
tion) . Recently, iGloveil ll2015ah identified the key chemical 
reactions for computing J„j,, and provided a reduced chem - 
ical network for this purpose. Subsequently. [Glove 3 (l2015tJ) 
investigated the impact of uncertainties in the chemical reac¬ 
tion rates, comparing their results with several chemical net¬ 
works employed in literature, including t he KROME pa ckage 
adopted here. In one-zone calculations, I Glove 3 (l2015lJl find 
that J„ii is at most a factor of a few higher than when using 
the KROME. Therefore, our results provide a lower limit on 
the estimates of Jem- However, we expect this difference to 
be further reduced in 3D simulations, as one-zone calcula¬ 
tions often underestimate the value of J„i, (iLatif et al.ll20lj . 
l2015h . 
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2.1 Sink Particles 


The simulations are evolved beyond the initial collapse 
by employing sink particles which represent the gravita¬ 
tionally bound objects. Such approach has been success¬ 
fully used both in smoothed particle hydrodynamics (SPHl 
and AMR code s (|Bate_et_alJ 199^ Krumholz et al.l l2004l : 


iFederrath et al.ll20ld : iLatif et all 20134 ). The general crite¬ 

ria used to create sink particles include assuming a density 
threshold, requiring the highest level of refinement, having 
converging flow, or considering the local minimum of grav¬ 
itational potential and/or a gravitationally bound region. 
In the literature, various combinations of these approaches 
are used and a detailed c ompa rison of these methods was 
made bv IFederrath et ^ (l2010l i. They found that a simple 
density threshold criterion to create sink particles is insuffi¬ 
cient in the regime of supersonic turbulence because o f loca l 
density enhancements. Recently, iBleuler fc Tevssi^ ll20l4l 
have presented a new method for sink particle creation in 
the Ramses code by employing a clump finding algorithm to 
more accurately locate the sites of sink formation by taking 
into account the tidal effects of surround ing gas distri b ution . 

We here employ the algorithm of IWang et al.l ll201(t) 
both for the creation of sinks, and their subsequent growth 
by accretion. Sinks are created only in the grid cells at 
the maximum refinement level, where, additionally, the Tru- 
elove criterion is violated, i.e., where the density exceeds the 
Jeans density, and in a region where the flow is convergent. 
Sinks are allowed to merge if they are formed within the 
accretion radius, i.e., the Jeans length in our case. Sinks 
are assigned masses so that a cell becomes Jeans stable 
after the subtraction of the sink mass. Sink velocities are 
computed based on m omen tum conservation. The wo rks of 
IFederrath et al.l (l201Cll ~) and iBleuler fc Tevssi^ (l2014l i sug¬ 
gest that this method of sink creation may create more 
sink particles. However, this may happen only in the regime 
where efficient fragmentation is expected, in the presence of 
supersonic turbulence or efficient cooling. We do not impose 
any accretion criterion such as the Bondi-Hoyle but allow 
the creation of additional sinks and let them merge within 
the accretion radius. We do not include any radiative feed¬ 
back. If one assumes that sinks m ay represent protostars , 
the stellar evolution calcul ations of ISchleicher et ^ ll2013li 
and iHosokawa et al.l ll2013li show that the ionising UV stel¬ 
lar feedback from the stars forming via rapid accretion, i.e., 
>0.1 Mo/yr remains very low even after they reach 10“^ Mo. 
Therefore, we expect that feedback is negligible as long as 
accretion rates are >0.1 Mo/yr. 


3 RESULTS 

We perform three dimensional cosmological simulations for 
three distinct massive primordial halos of a few 10^ Mo for 
various strengths of the BUV flux. The properties of the 
halos such as their masses, collapse redshifts, spins and the 
critical strength of UV flux are listed in table [T] In the fol¬ 
lowing, we discuss the structural properties of the halos and 
their temporal evolution. We also study the thermodynami¬ 
cal and physical properties of the halos. We distinguish here 
between the mass inflow rate, provided by the halo, and the 
accretion rate that occurs on a putative protostar. The mass 


Table 2. Eigenvalues of the inertia tensor for halo A. 


J 21 

ai 

a2 

as 

0 

0 

II 

7.9 

7.36 

4.37 

II 

0 

0 

3.0 

2.17 

2.14 

721 = 1000 

3.5 

3.4 

2.9 

721 = 40000 

7.7 

7.5 

6.7 


inflow rates are computed for the each intensity of the BUV 
flux and the mass accretion is followed by employing sink 
particles for a subset of halos. 


3.1 Structural properties of the halos 


According to our understanding of structure formation, the 
quantum fluctuations in the gravitational potential induced 
density perturbations. These perturbations collapsed via 
gravitational instabilities and merged with each other to 
form halos. Baryons fall in the DM potentials, get shock 
heated and this process continues until the virialization of 
the halo. We see that minihalos of a few times times 10^ Mq 
are formed at about z=20 in our simulations. Due to the 
presence of a BUV flux, these minihalos are unable to col¬ 
lapse via molecular hydrogen cooling. They continue to grow 
via mergers and accretion from the cosmic web until their 
virial temperature exceeds 10“* K. In figure [J] we show the 
structure of halo A at various cosmic times. This halo is 
formed at the intersection of filaments and is mainly fed 
by three main streams. The mass of the halo ~ 300 Myrs 
{z ~ 13) after the Big Bang is 1.4 x 10^ Mq and it undergoes 
a major merger with a halo of 5 x 10'’ Mq (i.e., mass ratio 
~1:2) about 100 Myrs later. The halo virializes at ~ z = 11, 
when its mass is well above the atomic cooling limit. In the 
meantime, it is continuously fed by gas flowing from the 
filaments and by minor mergers. 

Halo B goes through many minor mergers with mass 
ratios of about 1:10 and is also accreting gas via cosmic 
streams. Eventually, about 300 Myrs after the Big Bang, the 
halo becomes virialized, and has a mass of about 2.7x10^ Mq. 
This halo has a very different merger history from halo A 
as it has not experienced any major merger recently and is 
assembled ~130 Myr earlier than halo A. Similar behaviour 
is observed for all BUV flux strengths as the incident flux 
does not impact the large scale properties of the halos. Halo 
C experienced a major merger about 350 Myr after the Big 
Bang, similar to the history of halo A. The mass ratios of 
the merging halos is about 1:2 and the halo gets virialized 
about 50 million years later. 

Due to the presence of angular momentum in gas 
and halos, a central rotationally supported structure is ex- 
pected to form as a consequence of gravita tional collapse 
ilOh &: Haiman 20021: Volonteri fc Reej[2005h . However, as 
discussed by_[2b2LSi=i,!l2fi=^> triaxiality of dark mat¬ 


ter halos (| Allgood et aLlDOnfil i exerts gravitational torques 
and helps in transferring angular momentum via bars in 
bars instabilities, thus reducing rotational support compared 
to an axisymmetric case. In order to study the properties 
of the central regions, specifically their level of rotational 
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Figure 1. The ratio of enclosed DM mass to the enclosed baryonic mass is shown here for halo A, as a representative case. 


support, we have computed the ratio between the rota¬ 
tional velocity, Vroo and the velocity dispersion, cr. This ratio, 
'^roller, provides us the estimate of the disk scale height as 
H/R ~ cr/Vjot where H is the disk height, and R is the disk 
radius. Our estimates of Vroho- for halo A are shown in the 
bottom-right panels of figures |3l H] [5] and [6] for J2i=100, 500, 
1000 and 4 x 10**. Vroticr is ~ 6 in the central 10 pc of the 
halo for J 2 i = 100. Rotational support is higher at lower J 21 
values, as cooling in the core of the halo is efficient, because 
the fraction of molecular hydrogen (shown in the top-left 
panels) is higher. Vmilcr in fact becomes smaller for higher 
fluxes and approaches unity for the isothermal case. We es¬ 
timate that the aspect ratio (i.e. H/R) is ~ 0.25 for J 21 = 100 
and increases up to 0.5 for J 21 = 1000. The aspect ratio is 
almost unity for the isothermal case, i.e., for J 21 = 4 x lO'^. 
This suggests that a fat disk for all radiation field strengths 
is formed in the halo with aspect ratio higher than 0.1. Sim¬ 
ilar behaviour is observed for halo C (not shown). The value 
of Vror/cr is higher in the case of halo B (see figures, [71 [H] and 
m in general, but i t decreases as collap se approaches the 
isothermal state fcf. lOh fc Haimaal2002tl . It indicates that 
the scale height of the disk varies from halo to halo. 

We also calculated the components of the inertia tensor 
in the following way: 

h, = Yj nijiy] + z]) (1) 

/=1,A^ 

Ixy — — Y'^ ( 2 ) 

i=l.N 

The other components of the inertia tensor can be computed 


in a similar fashion. The diagonal elements of the inertia 
tensor are known as the moments of inertia. The eigenvalues 
ai, a 2 and aj of the inertia tensor for halo A are listed in 
tableland are computed at the same time at which Vroil o' 
is measured (i.e., last time evolution in figures[3l|4l[5]and[6]). 
For J 21 = 100, the eigenvalues of two components are almost 
equal but the value for the perpendicular component is 1.8 
times higher than the other two. The expected value for the 
perpendicular component of a thin disk is twice the values 
in the disk plane in agreement with our results. Similarly, 
the eigenvalues of two components are almost equal for J 21 = 
500 , 1000 and 40000 while the value of the perpendicular 
component is 1.4, 1.2 and 1.1 times the values in the disk 
plane respectively. This suggests that the disk is rather thin 
for I 21 = 100 and becomes thicker at higher fluxes. Similar 
trends are observed for the other two halos. 

We compute the rotational velocity by taking the ratio 
of the specific angular momentum i and the position vector 
r {Vrot ~ \{\I\A, the definition used in our previous analy¬ 
sis) and comp are it with the definition Ifl/lail given in the 
equation 13 of lReean fc Haehnelt 1 I 2 OO 9 II. These two defini - 
tions differ only by a factor of 1.3. Regan fc Haehn^ (l2009f l 
found similar differences, as expected in turbulent collapsing 
gas where rotation is not always well defined. We show the 
ratio of rotational to circular velocity in figure 1101 For the 
isothermal collapse case the ratio of Vrot/V^ir remains below 
unity except in the very centre (within the Jeans length) 
where we approach the resolution limit. The increase in Vjot 
above the V^irc by a factor of five for the weaker flux cases 
could be a consequence of various factors such as ellipsoidal 
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Figure 2. Average total (gas + DM) density for halo A illuminated by a BUV flux J 21 = 100. The columns show the structure of the 
halo at various cosmic times and the rows the total density averaged along the x, y &: z axis. The width of the density projections is 
given in physical units and the time in Myr after the Big Bang. In the middle column, the white circles highlight the main halo and the 
halo merging with it. 


rather than spherical geometry and large infall rates that 
keep the structure unsettled. This analysis confirms that 
rotational support is more significant in halos illuminated 
by low BUV fluxes compared with halos collapsing isother- 
mally. The higher rotational support in some halos can delay 
collapse and may even halt it if angular momentum is not 
transferred efficiently. In such a case, it may limit accretion 
onto the central object and have important implications for 
the final masses of stars. 


3.2 Thermodynamical properties 

The thermal properties of halo A for various strengths of 
BUV flux are shown in the figures [S] [d] [5] and [G] The tem¬ 
perature of the halo is a few thousand K when the halo 
mass is lower than the atomic cooling threshold. It increases 
through virialization and merger shocks until it reaches the 
atomic cooling limit. We note that in the presence of a BUV 
flux, irrespective of its strength, the temperature of the gas 
in the halo is ~ 8000 K after its virialization at scales larger 


than 10 pc. This is due to the lack of molecular hydrogen 
cooling as its abundance remains low, i.e. < 10“®. At smaller 
scales, H 2 molecules self-shield from the radiation depending 
on the strength of the ambient UV flux. For radiation field 
strengths < J 21 = 10“*, H 2 self-shielding becomes effective due 
to its (H 2 ) enhanced abundance and consequently the gas 
temperature decreases to a few hundred K within the central 
10 parsec. To further illustrate this effect, we show the H 2 
column density in figure [11] for various BUV fluxes for halo 
A. As expected, the H 2 column density becomes larger than 
10 '* for weaker fluxes. Consequently, H 2 self-shielding 
becomes effective and H 2 cooling occurs in the centre. On the 
other hand, for the strongest flux (i.e., J 21 = 4 x 10“*) the H 2 
column density remains lower than 10*^ cm“^, self-shielding 
remains ineffective and as a result the collapse is isothermal. 
We define the transition from atomic to molecular hydrogen 
cooling at H 2 abundance above 10“'*, which occurs between 
1-10 pc and is delayed by a few Myrs for the J 21 = 1000 case 
in comparison with J 21 = 100. The thermal evolution of halo 
B for J 21 = 100, 500 and 1000 is also very similar as shown 
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in figures [3 [8] and [9l Halo C (not shown) has also a very 
consistent behaviour. Their cores are cooled by molecular 
hydrogen with central temperatures of a few hundred K and 
have higher temperatures at larger radii. For the weakest 
flux cases, H2 cooling becomes important earlier than for 
the stronger flux cases. Variations are observed from halo to 
halo, e.g., in halo C the transition from atomic to molecu¬ 
lar state occurs at 0.4 pcil- S imilarly, variations from halo to 
halo are also observed fe.g.. lShang et aDl2010ll . 

In the case of the strongest flux, J21 = 4x 10"^, the forma¬ 
tion of molecular hydrogen remains suppressed throughout 
and the halo collapses isothermally with a central tempera¬ 
ture of 8000 K. The fraction of H2 in this case remains below 
10“^. Similar behaviour is observed for the other two halos 
for such high flux. Isothermal collapse occurs therefore only 
for significantly high values of J21, J21 > 10“*. The density in¬ 
side the halo increases during the collapse and the peak den¬ 
sity reached in the simulations is a few times 10“‘’ g/cm^'^. 
The density profile in halo A follows approximately an 
behaviour and becomes constant within the Jeans length. 
The small bumps in the profile indicate the presence of addi¬ 
tional substructure in the halo. The amount of substructure 
in the halos reduces with increasing flux strength. Similar 
features are observed for the other two halos. 


3.3 Mass inflow rates and masses at initial stages 

One of the prime objectives of the present work is to com¬ 
pute the mass inflow rates in the simulated halos for differ¬ 
ent thermodynamical conditions and assess whether they are 
sufficient to build a supermassive star. The latter requires 
mass accretion rates of a bout > 0.01 - 0.1 Mp,/vr to rapidly 
form a massive object (iBegelmarJ I 2 OI 0 I : iHosokawa et al.l 
l2013l : [Schleicher et al.ll2013l l. 

The temporal evolution of the mass inflow rates for halo 
A in the presence of various BUV flux strengths is shown in 
figures |3] (4] 0 and [6l The mass inflow rate is about 0.001 - 
0.01 Mg/yr when the halo mass is ~ 10® Mg. It increases up 
to ~ 0.1 Mg/yr when the halo mass is a few times lO’ Mg for 
all the strengths of BUV flux (J2i= 100, 500 and 1000). This 
is comparable to the theoretical free-fall rate, i.e. M ~ c^/G ~ 
0.1 Mg/yr (T/8000 where c, is the thermal sound speed. 
The inflow rate tends to be higher for J 21 = 1000 at all radii 
compared to the J 21 = 100 but differences are not significant. 
Similar values are found for J 21 = 500. In general, the small 
fluctuations in the inflow rates at various radii are due to the 
different density structure and the intermittent flows present 
in the halo (note that these are spherical averages). 

For halo B, the inflow rate for J21 = 100 is between 
0.01 - 0.1 Mg/yr while for J2i= 500 and 1000 it is about 
0.1 Mg/yr (see figures 0[8] and [9]). The drop in the inflow rate 
at 10 pc is due to the formation of an additional clump in the 
halo. Similar to halo A, the inflow rate for halo C for J2i = 
100, 500 and 1000 is ~ 0.1 Mg/yr. In all three halos and for all 
flux cases, the inflow rate declines within the Jeans length 
due to the higher thermal pressure. The variations in the 
inflow rates are correlated with the radial infall velocity of 

^ See lLatif et al.l ll2014h and lLatif et al.l ll2015l l where abundances 
of various chemical species are compared for the halos studied 
here. 


the halo. The typical infall velocities in these cases are 5-10 
km/s. The fluctuations in the infall velocity also indicate the 
onset of gravitational instabilities which may exert torques 
and help in the transfer of angular momentum. The infall 
velocity tends to increase with the strength of the BUV flux 
and becomes higher particularly for the isothermal collapse. 
Therefore, the mass inflow rate is larger for the isothermal 
collapse and is about 0.1-1 Mg/yr. 

We also estimated the total (DM -|- gas) enclosed mass 
in halo A. The enclosed mass within the central 10 pc of the 
halo is about a factor of a few higher for J 21 = 1000 compared 
to the J21 = 100 as shown in figures 0 [S] and 0 The mass 
profiles show a plateau around a few thousand solar masses 
which indicate the formation of small scale di sk, compara¬ 
ble to the estimates of Koi/o" and also found bv lReean et al.l 
ll2014h . This plateau becomes even more visible for halo B in 
the J 21 = 100 case. We also estimated the turbulent velocity 
for halo A in the presence of various strengths of UV flux 
and is shown in figures 000 and 0 The typical turbulent 
velocity is about 10 km/s and fluctuations in the turbulent 
velocity at various radii are related to the density structure 
and intermittent flows occurring. For the other two halos, 
we found similar trends both for the enclosed mass and tur¬ 
bulent velocities. Again, small variations in these quantities 
are observed from halo to halo. The gaseous structures (den¬ 
sity and temperature) within the central 10 pc of the three 
halos are shown in figures 112111311141 Either extended (halo 
A and halo B) or narrow funnels of dense gas connect to the 
halo centre, providing the high inflow rates we measured. 

3.4 Estimates of the masses of supermassive stars 

We selected two out of three halos, those with the highest 
inflow rates, to study the long term evolution of mass ac¬ 
cretion rates. We here only focus on the weaker BUV flux 
cases, where collapse is not fully isothermal. It has been 
shown that strong BUV fluxes lead to complete isothermal 
collapse where large accretion rates of about 1 Mg/yr can 
be obtained and massive objects of 10® Mg can be formed 
dRegan fc Haehn^ l2009l : iLatif et ^ l2013d : iBecerra et al.l 
l2015il . However. as explained earlier, these sites are very 
rare. As mentioned in the previous subsection, in halos A 
and C mass inflow rates of 0.1 Mg/yr can be obtained for 
J 21 > 500 and even for J 21 = 100. To assess if these rates can 
be maintained for long time scales, we employed sink parti¬ 
cles in our simulations and followed accretion onto them for 
about 30000 years after the formation of the first sink par¬ 
ticle. The sinks are formed at densities above 10^ and 
at the maximum refinement level. In general, no vigorous 
fragmentation is observed in either halo for all strengths of 
BUV flux. Only a single massive sink forms per halo. 

Figures [TS] and [13 show the time evolution of the gas 
density in the central 10 pc of halos A and C for J 21 = 1000 
after the formation of the sink particle. Note the highly 
anisotropic distribution of gas: the central sink is mainly 
fed by a narrow stream of gas, rather than by an extended, 
spherical or axisymmetric structure as in the case of a higher 
UV flux, where co llapse is isothermal (cf. figures 2 and 8 in 
iLatif et al.ll20133 l. 

The temporal evolution of the mass accretion rates onto 
the sinks for halo A, illuminated by various strengths of BUV 
flux, is shown in figure fTTl A mass accretion rate >0.1 Mg/yr 
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can be maintained for the simulated time for all cases. In 
the J 21 = 1000 case, the gas distribution at the time the sink 
particle is inserted is less dense and more extended, and it 
takes about 10000 yrs to move the smaller clumps into the 
centre and, after that, the inflow rate increases. Additionally, 
the small clumps that form, eventually merge in the centre. 
Therefore, in the last 5000 yrs of the simulated time, the 
mass accretion rate increases further for J21 = 1000 while the 
mass accretion rate for J 21 = 100 and 500 cases declines. 

The masses of the sinks after 30000 years, are 
4700, 6000 and 6770 Mo for J21 = 100, 500 and 1000, re¬ 
spectively. The mass of the sink is 2000 solar masses higher 
for J 21 = 1000 compared to the weakest flux case. In fact, the 
mass for J21 = 100 shows a plateau around 4700 solar masses 
and becomes almost constant for the last 10000 years. This 
is because of the higher rotational support as well as the 
presence of colder gas compared to the higher flux cases. 
The small plateau in the sink mass and a small decline in 
the mass accretion rate are also observed for J 21 = 500 in the 
last 5000 yrs. However, the mass of the sink seems to keep 
increasing for the J 21 = 1000 case. Therefore, we expect that 
the masses of the resultant star may be higher compared to 
the lower flux cases. 

We have also employed sink particles for halo C and fol¬ 
lowed the evolution for about 30000 yrs. The masses of the 
sinks and the mass accretion rates onto them are shown in 
figure [TS] The mass accretion rate starts to increase in the 
first 10000 yrs due to the sufficient mass supply in the central 
clump, but decline later. This behaviour is similar to halo A, 
but the decline is more significant. The sink mass starts to 
show a plateau earlier, at about 18000 yrs after the forma¬ 
tion of the sink. The mass accretion rates and the masses of 
the sink particles continue to increase for the rest of the two 
cases. The masses of the sinks are 4600, 6400 and 7000 Mq 
from the weaker to stronger fluxes, respectively (i.e. J21 = 
100, 500 and 1000). Similar to halo A, mass accretion rates 
remain above 0.1 Mg/yr for the simulated time, irrespective 
of the BUV field strength. The mass accretion rates show 
some decline for 721 = 100 case, but not for the J 2 i= 500 and 
1000 cases. Mass accretion rates may decline at later stages 
of evolution if rotational support is enhanced significantly, or 
efhcient fragmentation takes place inside the halo. However, 
for the duration of our simulation, we confirm >0.1 Mg/yr 
for the whole time. 

Overall, our simulations suggest that the mass accre¬ 
tion rates required to build a supermassive star can be ob¬ 
tained for UV fluxes below the critical value necessary for 
isothermal collapse. Particularly, the formation of supermas¬ 
sive stars seems plausible for J 21 = 1000 and 500, at least in 
one case. If accretion persists at the same rates, the masses 
of the sinks can reach lO'^ - 10 ^ solar masses. 


4 DISCUSSION AND CONCLUSIONS 

We have performed high resolution cosmological simulations 
for three distinct halos to study under which conditions the 
accretion rates necessary for the formation of supermassive 
stars in massive primordial halos of a few times 10 ^ Mq can 
be achieved. These supermassive stars are the potential cra¬ 
dles of direct collapse black holes which are one of the can¬ 
didates to explain the existence of high redshift quasars. 


The prime objective of this work is to explore the for¬ 
mation of supermassive stars in non-isothermal collapse. In¬ 
deed, the key parameter for supermassive st ar formation ap¬ 
pears to be an inflow r ate above 0.1 Mp/vr llBegelmanI 201 ^ 
jh^kawa et al.l I2OI3I : ISchleicher et al.l l2013l : iFerrara et al.l 


2 OI 4 II . Specificallv. we are interested in estimating whether 


the concept of a critical value, 7cm, of the BUV flux above 
which full isothermal collapse can be achieved through com¬ 
plete suppression of molecular hydrogen formation is neces¬ 
sary for black hole formation in the direct collapse scenario. 

We also study the structure of the inner regions of the 
halos to evaluate the role of rotational support. We pre¬ 
sume that halos are metal free and illuminated by moderate 
strengths of BUV flux. A pc-size fat disk is assembled in 
these halos as a consequence of gravitational collapse and 
the angular momentum of the infalling gas. The formation 
of such disk is irrespective of the merger history of the halo. 
The scale height of the disk varies from 0.25 - 1 depend¬ 
ing on the thermodynamical properties of the halo. Halos 
cooled by the molecular hydrogen have higher rotational 
support and an aspect ratio lower than the halos collaps¬ 
ing isothermally (lOh fc Haiman|l2002l : iLodato fc NataraianI 
I2OO6I : iBeeelmari et al.ll 20 o'^ 

Rotational support, however, does not, at least initially, 
hinder gas collapse in any of the simulated halos. In fact, our 
estimates for the mass inflow rates show that the halos irra¬ 
diated by moderate UV fluxes, J21 = 500 and 10 00 have typi¬ 
cal ma ss inflow rates of about 0.1 Mg/yr. Indeed. [Regan et al.l 
(l2014[l already noticed that sufficiently high accretion rates 
can be obtained for halos experiencing (anisotropic) fluxes 
of J21 ~ 1000. An inflow rate of ~ 0.1 Mg/yr can be achieved 
even for J 21 = 100 for halos A and C, while halo B has 
somewhat lower inflow rates. This may suggest that halos 
experiencing major mergers in their recent past may form 
supermassive stars more easily, for lower strengths of UV 
field, while perhaps this is not happeni ng in haloes f ormed 
via minor mergers and accretion. In fact. lMaver et all (l2014l ~) 
suggested that the merging of very massive galactic cores, 
> 10 ® Mo, may form a stable nuclear disk with very large gas 
inflow rates. They argue that merging of such systems leads 
to the enhanced inflow rates and helps in the formation of 
a massive black hole. Although the present work explores 
very different systems, the merging of a few times 10* Mq 
proto-galaxies, it hints that larger accretion rates seem to 
be obtained in the aftermath of a major merger, as mergers 
induce torques that remove angular momentum from gas, 
favoring inflows. However, our s mall sample does not allow 
us to draw strong conclusions. Ilnavoshi et al.l (l2015l l have 
proposed an alternative scenario which does not require the 
BUV to photo-dissociate molecular hydrogen but the colli¬ 
sions of protogalaxies with velocities of a few hundred km/s 
produce strong shocks with T > 10* K where subsequently 
gas collapses isobarically and H 2 gets collisionally dissoci¬ 
ated. 

To further assess whether such higher accretion rates 
can be maintained for long time scales, we have employed 
sink particles and followed accretion on them for 30000 years 
after the formation of the first sink particle for halos A and 
C. No strong fragmentation is found for these halos, irrespec¬ 
tive of the employed UV field, and a single star is formed per 
halo. A mass inflow ~ 0.1 Mg/yr can be maintained for long 
times, at least for the whole simulated time. We stopped the 
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simulations 30000 years after the formation of the first sink 
as evolving them further becomes computationally very ex¬ 
pensive. However, given these accretion rates, supermassive 
stars may reach 10“* - 10^ Mq within about one million years. 
Our results suggest that J 21 > 500 - 1000 may be sufficient 
to form supermassive stars, at least for some cases. The ob¬ 
served abundance of z > 6 quasars can be reproduced by 
the ab ove flux, extrapolating the findings of IPiikstra et ah] 
l|2014h . 

In our simulations, we employed sink particles to follow 
the mass accretion onto them, which introduce a character¬ 
istic length scale below which fragmentation is ignored. To 
assess fragmentation below this scale would require simu¬ 
lations resolving the collapse down to sub AU scales and 
it becomes computationally impractical to evolve them for 
longer times. This is due to the shorter time-steps involved in 
these calculations which cover a large range of spatial scales, 
from cosmological scales down to the scales of AU. However, 
the analysis of the density structure in these halos indicates 
that no vigorous fragmentation is expected, as H 2 cooling 
is only effective in the core of the halo and most of the gas 
outside the core is warm and no signs of substructures are 
present. Even if fr a gmen t ation occurs on very small scales, 
iLatif fc Schleichi^ 1120141 . [201^ show that clumps may mi¬ 
grate inward on short time scales. Moreover, viscous heating 
further dissociates molecular hydrogen, stabilises the disk 
and favours the formation of a massive central object. Based 
on the results from both th ese simulations and the comple - 
mentary analytical studies (iLatif fc Schleich^l2014l . l2015h . 
we argue that formation of supermassive stars in these con¬ 
ditions is possible. 

We also found that accretion onto the central star 
becomes highly anisotropic particularly for the case with 
J 21 = 1000 where the gas is fed to the central star by a sin¬ 
gle stream. Such behaviour seems to persist for 30000 yrs. 
This may have important implications for feedback by the 
central star or even by the black hole at later stages. The 
ionising radiation around a star/black hole would escape in 
the low density region while the gas supply to the central 
object may continue from a single dense stream even in the 
presence of feedback. 

The simulations presented in this work are based on 
a number of idealized assumptions such as halos are metal 
free and irradiated by an isotropic flux. In the future, sim¬ 
ulations self-consistently taking into account metal enrich¬ 
ment and radiation emitted by the nearby sources should be 
performed. Such simulations are mainly constrained by the 
exorbitant computational costs. 
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Figure 3. The thermodynamical and physical properties of halo A are shown for J 21 = 100. The top panels show the time evolution of 
gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, Vro?/cr, radial infall and 
turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Yrorl^^ =1- 
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Figure 4. The thermodynamical and physical properties of halo A are shown for J 21 = 500. The top panels show the time evolution of 
gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, Vro?/cr, radial infall and 
turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Vror/o" =1. 
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Figure 5. The thermodynamical and physical properties of the halo A are shown for J 21 = 1000. The top panels show the time evolution 
of gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, radial infall and 

turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Vror/o" =1. 
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Figure 6 . The thermodynamical and physical properties of the halo A are shown for J 21 = 4x 10*^. The top panels show the time evolution 
of gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, Vror/cr, radial infall and 
turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Vror/o" =1. 
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Figure 7. The thermodynamical and physical properties of the halo B are shown for J 21 = 100. The top panels show the time evolution 
of gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, Vror/o"? radial infall and 
turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Vrorl^^ =1- 
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Figure 8 . The thermodynamical and physical properties of the halo B are shown for J 21 = 500. The top panels show the time evolution 
of gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, radial infall and 

turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Vror/o" =1. 
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Figure 9. The thermodynamical and physical properties of the halo B are shown for J 21 = 1000. The top panels show the time evolution 
of gas density, temperature, H 2 fraction and of the total enclosed mass. The time evolution of the inflow rates, radial infall and 

turbulent velocities is shown in the bottom panels. The different line styles represent various cosmic times. The profiles are spherically 
averaged and radially binned. The vertical dashed black line shows the Jeans length and the horizontal dashed line represents Vror/o" =1. 
























Black hole formation in the early universe 



Radius [pc] 


Figure 10. The ratio of rotational to circular velocity is shown here for various BUV fluxes for halo A, as a representative case. 
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Figure 11. The H 2 column density for various BUV fluxes is shown here for halo A, as a reference. 
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Figure 12. Average density (top panels) and temperature (bottom panels) weighted by the gas density along the y-axis for Halo A. 
They are plotted for the BUV of strengths J 21 =100, 500 &; 1000 and show the central 10 pc region of the halo when the simulation 
reaches the maximum refinement level prior to the formation of a sink particle. 
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Figure 13. Average density (top panels) and temperature (bottom panels) weighted by the gas density along the y-axis for Halo B. 
They are plotted for the BUV of strengths J 21 =100, 500 Sz 1000 and show the central 10 pc region of the halo when the simulation 
reaches the maximum refinement level prior to the formation of a sink particle. 
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Figure 14. Average density (top panels) and temperature (bottom panels) weighted by the gas density along the y-axis for Halo C. 
They are plotted for the BUV of strengths J 21 =100, 500 &; 1000 and show the central 10 pc region of the halo when the simulation 
reaches the maximum refinement level prior to the formation of a sink particle. 
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Figure 15. The time evolution of the density averaged along the y-axis is shown in the central 10 pc of the halo A for J2i=1000 after 
the formation of the sink particle. The time after the formation of the sink particle in each panel is shown in years and the white dot 
depicts the position of the sink particle. 
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Figure 16. The time evolution of the density averages along the y-axis is shown in the central 10 pc of the halo C for J2i=1000 after 
the formation of the sink particle. The time after the formation of the sink particle in each panel is shown in years and the white dot 
depicts the position of the sink particle. 
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Figure 17. The masses of the sink particles and their mass ac¬ 
cretion rate are shown against time for halo A. The strength of 
incident field is shown by different colors and line styles as de¬ 
scribed in the legend. 
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